/* pde.h */
#ifndef EQHEAT_H_
#define EQHEAT_H_

#include <armadillo>
#include <vector>

typedef double (*ptr2func1arg)(double);				// real function of one variabe
typedef double (*ptr2func2arg)(double, double);		// real function of two variale
typedef ptr2func1arg f1v;							// short name for real function of one variabe
typedef ptr2func2arg f2v;							// short name for real function of two variabe

class EQHeat
{
private:
	double BC_L;   				// container for BC_left
	double BC_R;				// container for BC_right
	std::vector<arma::vec> u;	// container for discrete values of u(x,t)
	double a;					// diffusion coefficient
	double xlimits[2];			// space domain		
	double tlimits[2];			// time domain
	double hx;					// step in space domain
	double ht;					// step in time domain
	unsigned int Nx;			// number of points (length) in space domain
	unsigned int Nt;			// number of points (length) in time domian


	/* the Heat equation is u_{x} + \frac{1}{a^2}u_{xx} = 0, 0 <= x <= l, t>0; 
	 * with boundary conditions ans initial condition:
	 * u(0,t) = g_1(t) -------------> BC_left
	 * u(l,t) = g_2(t) -------------> BC_right
	 * u(x,0) = g_3(t) -------------> IC
	 */
	 
	f1v _BC_left;
	f1v _BC_right;
	f1v _IC;

	/* By the implict difference method, it can obtain an algebraic equation 
	 *
	 * 							Au^{n+1}=u^n
	 *
	 * where A is a sysmtric banded matrix,
	 *
	 * 	   [ 1+2d     -d                       ]      
	 *     [   -d   1+2d     -d                ]		
	 * A = [    	 ...    ...    ...         ] /a/a= toeplitz([1+2d -d ... ])
	 * 	   [                 -d   1+2d      -d ]		
	 *     [                        -d    1+2d ]		
	 *
	 * d = a*a*ht/hx/hx;
	 */
	double d;
	arma::mat A; 

	unsigned int compute_Nx();
	unsigned int compute_Nt();
	void init();

public:

	/* parameter settors */
	void set_a(double);
	void set_xlimits(double lower, double upper);
	void set_tlimits(double lower, double upper);
	void set_BC_left(f1v);
	void set_BC_right(f1v);
	void set_IC(f1v);
	void set_hx(double);
	void set_ht(double);
	void solve();

	void save_data();
	void save_data(const char* filename) const;
	void draw_data();
	void draw_data(const char* filename, std::vector<double> time) const;
};

#endif